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Abstract 

This  work  deals  with  the  control  of  polymeric  fuel  cells.  It  includes  a  linear  analysis  of  the  system  at  different  operating  points,  the 
comparison  and  selection  of  different  control  structures,  and  the  validation  of  the  controlled  system  by  simulation.  The  work  is  based  on  a 
complex  non  linear  model  which  has  been  linearised  at  several  operating  points.  The  linear  analysis  tools  used  are  the  Morari  resiliency  index, 
the  condition  number,  and  the  relative  gain  array.  These  techniques  are  employed  to  compare  the  controllability  of  the  system  with  different 
control  structures  and  at  different  operating  conditions.  According  to  the  results,  the  most  promising  control  structures  are  selected  and  their 
performance  with  PI  based  diagonal  controllers  is  evaluated  through  simulations  with  the  complete  non  linear  model.  The  range  of  operability 
of  the  examined  control  structures  is  compared.  Conclusions  indicate  good  performance  of  several  diagonal  linear  controllers.  However,  very 
few  have  a  wide  operability  range. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Proton  exchange  membrane  fuel  cells  (PEMFC),  also 
known  as  polymer  electrolyte  fuel  cells,  are  being  actively 
developed  for  use  in  cars,  buses,  as  well  as  for  a  very  wide 
range  of  portable  applications,  and  also  for  combined  heat 
and  power  systems  [1,2].  They  are  regarded  as  ideally  suited 
for  transportation  applications  due  to  its  high  power  density, 
compactness,  lightweight,  low-operating  temperature  which 
permits  fast  start-up,  solid  electrolyte,  long  cell  and  stack 
live,  low  corrosion,  and  higher  efficiencies  compared  to  heat 
engines  [3,4]. 

A  lot  of  works  are  dedicated  to  model  polymeric  fuel 
cells,  as  summarised  by  Yao  et  al.  [5].  However,  most  of 
these  models  are  static  models.  Some  authors  introduce 
electrochemical  time  constants  in  their  models  [6],  but  they 
do  not  model  the  flow  dynamics  of  the  whole  system.  With 
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the  aim  to  study  the  flow  dynamics  of  polymeric  fuel  cells, 
Pukrushpan  et  al.  [7]  presented  a  control  oriented  model 
which  includes  the  transient  phenomena  of  the  compressor, 
the  manifold  filling  dynamics  (both  anode  and  cathode), 
reactant  partial  pressures,  and  membrane  humidity.  This 
model  proposed  by  Pukrushpan  et  al.  has  been  the  base  for 
the  model  used  in  this  work  although  some  modifications  are 
introduced. 

In  1998,  Yang  et  al.  [8]  described  the  control  challenges  in 
fuel  cell  vehicle  development.  After  that,  some  studies  have 
addressed  these  challenges  related  to  the  control  of  polymeric 
fuel  cells  [4,9].  This  literature  has  been  considered  to  define 
the  multiple  objectives  of  the  control  problem  addressed  in 
this  work. 

PI  based  controllers  are  proposed  by  some  authors  but  a 
complete  study  of  the  multiple  input  multiple  output  (MIMO) 
control  problem  is  not  done.  In  this  work,  different  control 
structures  are  analysed  and  compared  taking  into  account 
the  interactions  between  loops  and  the  directionality  of  the 
system. 
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Nomenclature 

Afc 

active  area  per  cell 

Cy 

water  concentration 

ht 

stack  load  current 

i 

current  density 

ASm,out 

fluidynamic  constant  at  supply  manifold  out¬ 
put 

ACa,out 

fluidynamic  constant  at  the  cathode  output 

Ann,  out 

fluidynamic  constant  at  the  return  manifold 
output 

Aan.out 

fluidynamic  constant  at  the  anode  output 

m\ 

mass  of  component  i 

n 

number  of  fuel  cells 

Pea. 

cathode  pressure 

Po2 

cathode  oxygen  partial  pressure 

m2 

anode  hydrogen  partial  pressure 

Psat 

saturation  pressure  of  water  at  Tfc 

A p 

anode-cathode  pressure  difference 

Tic 

fuel  cell  temperature 

fin 

membrane  thickness 

Vfc 

unique  fuel  cell  voltage 

fim 

the  compressor  voltage 

^st 

stack  voltage 

mass  flowrate  of  component  i 

•*442,  an,  in 

the  mass  fraction  of  the  hydrogen 

3;02,ca,in 

the  molar  fraction  of  input  air 

Greek  letters 

0ca,in 

the  relative  humidity  of  the  input  air 

Xq2 

oxygen  utilisation 

fim 

membrane  water  content 

To  use  fuel  cells  in  specific  applications,  electric  con¬ 
ditioning  systems  are  normally  added.  The  complexity  and 
performance  of  these  power  electronic  systems  vary  and 
consequently  the  stress  suffered  by  the  fuel  cell  varies  too. 
Most  literature  works  analyse  the  ability  of  the  fuel  cell 
control  system  to  follow  changes  in  the  current  load,  and 
assume  large  current  excursions.  However,  thanks  to  the  con¬ 
ditioning  systems,  these  excursions  can  be  minimised,  and 
the  fuel  cell  maintained  close  to  a  nominal  operating  point. 
In  this  small  region  behaviour  is  closer  to  linear  behaviour.  In 
automotive  applications,  for  example,  fuel  cells  can  be  used 
in  hybrid  configurations  with  batteries  allowing  the  cells  to 
be  operated  at  steady  state  with  high-power  output  under 
optimised  conditions,  and  with  the  peak  power  demands 
being  met  by  batteries  [10].  Alternatively,  the  fuel  cells  can 
be  used  as  the  sole  power  supply  and  sudden  load  application 
are  a  possible  scenario.  For  this  reason,  in  this  work,  large 
changes  in  the  current  load  as  well  as  small  excursions  from 
a  nominal  operating  point  are  taken  into  account. 

The  rest  of  the  paper  is  structured  as  follows:  in  Section 
2,  the  work  methodology  is  described,  Section  3  contains  the 


description  of  the  model,  in  Section  4  the  control  objectives 
are  explained,  in  Section  5  the  linear  analysis  is  done  and  the 
preferred  control  structures  are  identified,  in  Section  6  the 
preferred  control  structures  are  validated  through  simulation, 
in  Section  7  the  best  control  possibilities  are  discussed,  and 
finally,  in  Section  8  the  main  conclusions  are  explained.  Two 
appendices  recall  some  of  the  model  equations  and  in  the 
Nomenclature  Section  the  different  variables  that  appear  in 
the  work  are  defined. 


2.  Methodology 

This  work  is  based  on  a  non  linear  model  of  a  PEMFC  pro¬ 
posed  by  Pukrushpan  et  al.  [7],  which  has  been  linearised  at 
different  operating  points.  MIMO  linear  models  are  obtained 
in  the  form  of  state  space  matrices  A,  B,  C,  and  D. 

The  work  faces  the  control  of  two  output  variables  us¬ 
ing  two  of  the  free  input  variables.  Therefore,  2x2  control 
structures  are  considered.  In  this  work,  a  control  structure  is 
understood  as  the  set  of  inputs  chosen  to  control  the  system, 
as  well  as  the  input-output  pairing.  Since  there  exist  more 
than  two  possible  manipulated  variables,  there  exist  different 
pairs  of  input  variables  to  control  the  system.  From  a  linear 
system  that  includes  all  the  input  variables,  the  different  2x2 
linear  systems  corresponding  to  the  different  input  pairs  are 
derived  removing,  in  each  case,  the  columns  of  the  rejected 
inputs. 

MIMO  linear  systems  can  be  analysed  using  different 
analysis  tools.  These  tools  are  mathematic  operators  applied 
to  the  squared  transfer  functions  of  the  linear  system  that 
give  relevant  information  such  as  stability,  controllability, 
sensitivity,  robustness,  etc.  To  design  the  control  of  MIMO 
systems,  one  of  the  most  important  tasks  is  the  selection 
of  the  set  of  manipulated  variables  to  control  the  system. 
Different  controllability  indexes  can  help  the  designer 
choose  this  variables  set. 

In  the  literature,  three  different  controllability  indexes  are 
frequently  used  [11]:  the  Morari  resiliency  index  (MRI),  the 
condition  number  (CN),  and  the  relative  gain  array  (RGA). 
These  indexes  are  used  in  this  work  to  study  the  interaction 
between  control  loops  and  the  sensitivity  of  the  controlled 
system. 

The  singular  value  decomposition  is  a  numerical  al¬ 
gorithm  useful  in  analysing  the  multivariable  aspect  of 
the  gain  matrix,  giving  the  input  and  output  directions  for 
which  gains  are  maximum  and  minimum.  The  MRI  is  the 
smallest  singular  value  of  the  open-loop  transfer  function. 
It  is  the  poorer  gain  of  the  process,  poorer  sensitivity,  which 
corresponds  to  specific  input  and  output  directions.  MRI  is 
one  of  the  controllability  indexes  calculated  and  analysed 
in  this  work.  Control  structures  with  large  MRI  over  the 
frequency  range  of  interest  are  preferred. 

The  second  controllability  index  analysed  is  the  CN.  It  also 
comes  from  the  singular  value  decomposition  of  the  transfer 
function.  The  CN  is  the  ratio  of  the  maximum  singular  value 
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to  the  minimum  singular  value  and  it  is  typically  used  for 
the  control  structure  selection.  It  provides  a  numerical  indi¬ 
cation  of  the  sensitivity  balance  in  a  multivariable  system. 
Large  condition  numbers  indicate  unbalanced  sensitivity  and 
also  sensitivity  to  changes  in  process  parameters.  Therefore, 
structures  with  small  CNs  are  preferred. 

The  third  index  analysed  is  the  RGA.  RGA  is  used  to 
determine  the  interaction  among  control  loops  in  a  multivari¬ 
able  process.  It  is  defined  as  the  ratio  of  the  open-loop  gain 
for  a  selected  output  when  all  the  other  loops  of  the  process 
are  open,  to  its  open-loop  gain  when  all  the  other  loops  are 
closed.  RGA  of  a  complex  non-singular  matrix  M  is  calcu¬ 
lated  as  indicated  in  Eq.  (1),  where  X  denotes  element  by 
element  multiplication  (Hadamard  product). 

RGA(M)  =  MX(M“ 1  )T  ( 1 ) 

Pairings  that  have  RGA  close  to  unity  matrix  at  frequen¬ 
cies  around  bandwidth  are  preferred.  This  rule  favours  min¬ 
imal  interaction  between  loops,  which  means  independence 
of  the  loops.  Being  the  loops  independent,  stability  problems 
caused  by  interaction  are  prevented.  Numbers  around  0.5  in¬ 
dicate  relevant  interaction.  The  RGA  indicates  other  useful 
control  properties  [11].  One  of  the  most  important,  structures 
with  large  RGA  elements  around  the  bandwidth  frequency  are 
difficult  to  control  because  of  sensitivity  to  input  uncertainty. 

In  this  work,  these  three  indexes  are  used  to  characterise 
the  controllability  of  polymeric  fuel  cells.  They  are  used  to 
compare  different  control  structures  for  a  PEMFC  and  a  con¬ 
trol  structure  at  different  operating  conditions. 

Normally,  only  the  steady  state  value  of  these  controlla¬ 
bility  indexes  is  regarded.  However,  their  analysis  in  the  fre¬ 
quency  domain  is  important  [11].  Consequently,  in  this  work, 
the  three  indexes  are  analysed  in  a  wide  frequency  range. 
The  frequency  range  of  interest  is  given  by  the  bandwidth 
frequency.  The  bandwidth  frequency  is  normally  defined  as 
the  frequency  up  to  which  control  is  effective  and  it  depends 
on  the  controller. 

The  main  goal  of  this  work  is  the  linear  analysis  of  the  sys¬ 
tem.  However,  the  PEMFC  is  not  linear.  Linearised  model 
behaviour  cannot  be  extrapolated  when  the  system  makes 
large  excursions  from  the  nominal  operating  point.  To  prove 
the  capability  of  the  preferred  control  structures  to  drive  the 
system  correctly,  simulations  based  on  a  non  linear  system 
are  necessary.  Because  of  that,  after  the  controllability  anal¬ 
ysis  in  Section  5,  an  analysis  based  on  the  non  linear  model 
simulations  is  done  (Section  6). 


3.  Model 

3.1.  Nonlinear  model 

Results  of  this  work  are  based  on  a  non  linear  model  pro¬ 
posed  by  Pukrushpan  et  al.  [7].  The  model  was  developed 
specifically  for  control  with  the  idea  of  avoiding  unnecessary 
detail.  On  the  other  hand,  special  attention  is  given  to  the  air 


circuit  subsystem  and  a  quite  accurate  compressor  model  is 
included. 

Important  simplifications  are  done  because  for  control  pur¬ 
poses,  excessive  detail  is  not  necessary.  Because  of  that,  spa¬ 
tial  variations  are  not  included  and  constant  properties  are 
assumed  in  all  volumes.  Only  time  derivatives  are  present. 
With  respect  to  the  considered  dynamics,  the  model  neglects 
the  fast  dynamics  of  electrochemical  reactions  and  electrode 
electricity.  Temperature  is  treated  as  a  constant  parameter 
because  its  slow  behaviour  (time  constant  of  about  102  s) 
permits  that  it  is  regulated  with  its  own  (slower)  controller. 

The  fuel  stack  model  contains  four  interacting  parts:  the 
stack  voltage,  membrane  hydration,  the  cathode,  and  the 
anode. 

Stack  voltage,  nst,  is  calculated  as  a  function  of  stack  cur¬ 
rent,  /st,  cathode  pressure,  pCSL,  reactant  partial  pressures,  po2 
and  pw2,  fuel  cell  temperature,  7fc,  and  membrane  humidity. 
Identical  behaviour  of  each  cell  is  assumed  and  stack  voltage 
is  calculated  as  the  individual  cell  voltage,  Vfc ,  per  the  num¬ 
ber  of  cells,  n.  The  cell  voltage  has  four  terms,  as  can  be  seen 
in  the  following  equation 

^fc  —  E  —  nact  ^ohm  t’conc  (2) 


where  E ,  the  open  circuit  voltage,  is  a  function  of  the  fuel  cell 
temperature  and  hydrogen  and  oxygen  partial  pressures;  nact, 
the  activation  overvoltage,  depends  on  the  current,  the  tem¬ 
perature,  and  the  oxygen  partial  pressure;  n0hm>  the  ohmic 
overvoltage,  is  proportional  to  the  stack  current  and  has  a 
proportionality  constant  strongly  dependent  on  membrane 
humidity  and  temperature;  and  vconc,  the  concentration  over¬ 
voltage,  depends  on  the  current  and  oxygen  partial  pressure. 
In  Appendix  A,  detailed  expressions  are  shown. 

Membrane  hydration  captures  the  effect  of  water  transport 
across  the  membrane.  Both  water  content  and  mass  flow  are 
assumed  to  be  uniform  over  the  surface  area  of  the  membrane, 
and  are  functions  of  stack  current  and  relative  humidity  of 
the  gas  in  the  anode  and  cathode.  Detailed  hydration  model 
is  given  in  Appendix  B. 

Inside  the  cathode  volume,  the  mass  flow  continuity  is  used 
to  balance  the  mass  of  the  three  elements  (oxygen,  nitrogen, 
and  water).  Some  assumptions  are  made:  ideal  gases,  cathode 
temperature  equal  to  the  stack  temperature,  same  properties 
of  the  gas  exiting  the  cathode  and  the  gas  inside  the  cathode, 
and  flow  channels  and  backing  layer  lumped  into  one  volume. 
Equations  for  the  mass  time  derivatives  are: 


W02  ,out  Wo2,  reacted 


dmN2 
d  t 


—  WN2,in  WN2,out 


(3) 

(4) 


dm 


w,ca 


d  t 


—  fkv,ca,in  Wv,ca,out  +  Wy  ,ca,gen  T  Wv,membr 


n  l§t 

Wo2,  reacted  =  T/q2  x 
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fkv,ca,gen  —  47v  x 


nht 

~2F 


where  W-m  and  Wout  are  the  flow  rates  entering  and  exiting  the 
cathode,  Wq2, reacted  is  the  mass  of  oxygen  reacted  per  unit  of 
time,  Wv,ca,gen  is  the  mass  of  water  generated  by  the  reaction 
per  unit  of  time,  Wv,membr  is  the  mass  of  water  that  crosses 
the  membrane  per  unit  of  time,  Mo2  and  Mv  the  molar  masses 
of  oxygen  and  water,  and  F  the  Faraday  constant. 

The  air  entering  the  cathode  is  impelled  by  a  compressor. 
The  modelled  system  includes  a  compressor  with  a  dynamic 
part  and  a  static  part  the  behaviour  of  which  is  read  from  an 
experimental  compressor  map.  Supply  and  return  manifolds, 
cooler  and  humidifier  are  also  included.  Supply  and  return 
manifolds  are  lumped  and  modelled  as  unique  volumes.  The 
supply  manifold  pressure  is  governed  by  mass  continuity  and 
energy  conservation  equations  assuming  ideal  gas.  The  re¬ 
turn  manifold  pressure  is  governed  by  the  mass  continuity 
and  the  ideal  gas  law  through  isothermal  assumptions.  An 
ideal  cooler  maintains  the  temperature  of  the  air  entering  the 
stack  at  a  fixed  point.  The  pressure  drop  across  the  cooler  is 
negligible.  It  is  assumed  that  a  static  humidifier  gives  to  the 
air  the  desired  relative  humidity  before  entering  the  stack  and 
the  amount  of  water  injected  is  calculated. 

Flow  rates  from  one  volume  to  another  are  calculated 
in  function  of  the  upstream  and  downstream  pressures.  The 
same  linear  behaviour  is  assumed  at  the  supply  manifold  out¬ 
put,  cathode  output,  and  return  manifold  output.  In  all  cases, 
the  following  expression  applies,  where  ^upstream  is  the  up¬ 
stream  pressure  and  /^downstream  is  the  downstream  pressure: 


W  —  K nozzle  x  (/^upstream  /^downstream)- 


At  the  anode  side,  entering  hydrogen  comes  from  a  pres¬ 
surised  tank.  Since  pressurised  hydrogen  is  used,  the  hydro¬ 
gen  flow  is  assumed  to  be  a  manipulated  input  variable.  Al¬ 
though  the  model  used  in  this  work  is  based  on  the  model 
presented  by  Pukrushpan  et  al.  [7],  an  important  modifica¬ 
tion  is  done:  the  addition  of  an  exit  flow  in  the  anode.  This  is 
necessary  to  keep  the  system  flooded  with  gas  and  to  improve 
the  power  demand  transient  responses  [8] .  To  model  this  flow, 
a  linear  behaviour  as  indicated  in  Eq.  (8)  is  assumed. 

The  anode  flow  model  is  quite  similar  to  the  cathode  flow 
model.  Hydrogen  partial  pressure  and  anode  flow  humidity 
are  determined  by  balancing  the  mass  of  hydrogen  and  water 
in  the  anode. 


d 

dt 
dm 


=  Wu2,m  - 


w,an 


dt 


=  W, 


v,an,m 


Wn2  ,out  IFh2, reacted 

—  fFv,an,out  fFy^membr 


(9) 

(10) 


where  W-m  and  Wout  are  the  flow  rates  entering  and  exiting 
the  anode  and  Wu2, reacted  is  the  mass  of  hydrogen  reacted  per 
unit  of  time. 

All  together,  it  results  in  a  model  of  nine  states:  the  angular 
velocity  of  the  compressor,  the  mass  of  oxygen  in  the  cathode, 
the  mass  of  nitrogen  in  the  cathode,  the  mass  of  water  vapour 


in  the  cathode  (saturated  for  some  operating  conditions),  the 
mass  of  hydrogen  in  the  anode,  the  mass  of  water  vapour 
in  the  anode,  the  pressure  in  the  supply  manifold,  the  mass 
of  air  in  the  supply  manifold,  and  the  pressure  in  the  return 
manifold.  In  case  of  water  saturation,  eight  states  are  contem¬ 
plated  instead  of  nine  because  the  mass  of  water  vapour  in 
the  cathode  is  removed.  More  detail  of  the  non  linear  model 
can  be  found  in  [7]. 

3.2.  Model  linearisation 

The  non  linear  model  described  above  has  been  linearised 
at  different  operating  points.  SIMULINK®  linearisation  tools 
have  been  used  to  obtain  the  state  space  matrices  of  the  sys¬ 
tem. 

MRI  and  CN  are  scale  dependent.  Because  of  that,  the  lin¬ 
ear  models  have  been  scaled.  One  of  the  controlled  outputs 
is  the  difference  of  pressure  between  anode  and  cathode,  A p, 
with  a  maximum  accepted  variation  of  0. 1  bar.  This  has  been 
used  as  the  scaling  parameter.  For  the  rest  of  input  and  output 
variables,  scaling  has  been  done  assuming  a  maximum  vari¬ 
ation  of  10%.  Hence,  the  scaled  variables  are  the  non  scaled 
increments  divided  by  the  maximum  increments. 

3.3.  Operating  conditions 

The  modelled  fuel  cell  has  a  number  of  cells  ft  =  381,  an 
active  area  per  cell  Afc  =  280  cm2,  and  an  operating  temper¬ 
ature  7fc  =  80  °C.  Linearisation  of  the  system  has  been  done 
at  four  different  operating  points  that  have  been  called  OP1, 
OP2,  OP3,  and  OP4.  In  all  cases,  the  parameters  of  the  model 
predict  full  humidified  cathode  and  the  water  vapour  in  the 
cathode  is  saturated  [7]. 

For  OP1,  a  net  power  of  Pnet  =  37,400  W  is  required.  The 
net  power  is  the  power  given  by  the  fuel  cell  minus  the  power 
consumed  by  the  compressor.  The  desired  amount  of  power 
can  be  obtained  at  different  load  currents.  However,  there 
is  a  minimum  load  current  for  which  this  power  can  be  ob¬ 
tained,  and  its  value  is  7st  =  175  A.  With  the  aim  of  having 
OP1  with  the  minimum  hydrogen  consumption  (minimum 
hydrogen  reaction  if  hydrogen  is  recycled),  an  operating  point 
with  7st  =  175  A  is  chosen.  OP1  is  then  defined  by  7st  =  175  A, 
vcm  =  158  V,  Wan?in  =  1.0134  g  s_1,  giving  an  output  stack 
voltage  of  i>st  =  242.75  V,  pCSL  =  1.99  bar,  A p  =  —0.0068  bar, 
and  Pnet  =  37,400  W.  With  these  conditions  the  system  is  op¬ 
erating  in  the  linear  zone  of  the  polarisation  curve,  with  a 
current  density  i  =  0.62  A  cm-2. 

OP2  and  OP3  have  the  same  output  voltage  than  OP1 
and  the  same  A p,  but  different  current  loads.  OP2  is  de¬ 
fined  by  7st  =  150  A  (closer  to  the  activation  zone),  and  has 
vcm  =  137  V,  Wan5in  =  0.86gs_1,  and  pCa  =  1-82  bar;  OP3  is 
defined  by  7st  =  200A  (closer  to  the  concentration  zone), 
and  has  vcm  =  183  V,  Wan?in  =  1 . 1 8  g  s_  1 ,  and  pCSL  =  2.21  bar. 
Therefore,  OP2  and  OP3  are  operating  points  the  controlled 
system  will  pass  through  with  load  currents  of  150  and  200  A. 

Finally,  OP4  is  defined  fixing  the  net  power  at  37,400  W 
and  7st  =  200A.  The  corresponding  compressor  voltage  is 
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vcm  =  1 30 . 24  V,  vst  =  20 1 . 85  V,  and  j?ca  =1.73  bar.  OP4  has 
the  same  net  power  as  OP1  but  a  higher  hydrogen  consump¬ 
tion.  It  has  been  chosen  to  compare  the  controllability  of  the 
system  when  it  is  operated  at  optimal  and  non-optimal  con¬ 
ditions  with  respect  to  the  hydrogen  consumption. 

4.  Analysis  of  the  control  problem  and  control 
objectives 

The  first  control  objective  addressed  in  this  work  is  to 
maintain  vst  to  the  setpoint  value.  On  the  other  hand,  in  or¬ 
der  to  prevent  membrane  damage,  the  difference  of  pressure 
between  anode  and  cathode  A p  has  to  be  small  [8].  The  ad¬ 
equate  A p  value  will  depend  on  the  membrane  support  and 
on  the  age  of  the  fuel  cell.  However,  since  minimum  A p  will 
favour  the  membrane  life  time,  the  second  control  objective 
addressed  is  to  maintain  A p  close  to  zero. 

The  inputs  of  the  system  are  the  oxygen  molar  fraction 
of  the  dry  input  air,  yo2,ca,in>  the  compressor  voltage,  ncm, 
the  flow  of  hydrogen  Wan, in,  the  relative  humidity  of  the  in¬ 
put  air,  0ca,in>  the  mass  fraction  of  the  hydrogen  (hydrogen 
versus  water)  VH2,an,in>  and  the  fluidynamic  constants  at  the 
supply  manifold  output,  cathode  output,  return  manifold  out¬ 
put,  and  anode  output,  Ksmoui,  KCSi0Ui,  A7m,out>  and  an, out? 
respectively.  Some  of  these  inputs  may  not  be  manipulated 
variables  in  real  applications,  as  will  be  discussed  in  Section 
7.  However,  their  inclusion  in  the  controllability  analysis  may 
give  important  information  about  their  influence. 

As  has  been  explained,  there  are  two  output  variables  to 
be  controlled  and  nine  different  inputs.  In  consequence,  there 
are  multiple  control  structures  which  are  possible.  Character¬ 
isation  and  comparison  between  them  is  done  in  Sections  5 
and  6. 

In  the  literature,  the  importance  of  the  control  of  oxygen 
excess  ratio  Xq2  stressed  [4].  Abrupt  changes  of  this  vari¬ 
able  should  be  avoided  in  order  to  prevent  abrupt  changes 
in  the  oxygen  partial  pressure  and  oxygen  starvation.  The 
control  of  Xo2  instead  of  the  control  of  vst  is  also  studied  in 
Section  6. 

From  this  point  on,  control  structures  are  designed  by  a 
pair  of  control  variables,  and  it  is  understood  that  the  first  of 
them  controls  nst  and  the  second  of  them  controls  A p. 

5.  Results  of  the  linear  analysis 

The  results  of  the  analysis  of  MRI,  CN  and  RGA  in  the  fre¬ 
quency  range  between  10-2  and  102  rad  s_1  are  summarised 
in  this  section.  Frequencies  larger  than  102rads_1  are  not 
analysed  because  it  is  considered  that  it  is  not  convenient  for 
the  system  to  follow  such  rapid  frequencies  and,  therefore, 
bandwidth  frequency  will  not  be  larger  than  102  rad  s_1 . 

Since  two  manipulated  variables  have  to  be  selected  from 
nine  possibilities,  there  are  36  possible  sets,  or  control  struc¬ 
tures.  All  36  sets  have  been  compared  using  the  three  indexes 


Frequency  (rad/s) 


Fig.  1.  RGA  at  operating  point  one. 


Fig.  2.  RGA  at  operating  point  four. 


in  the  frequency  domain.  However,  in  Figs.  1-6,  only  the  in¬ 
dexes  with  best  behaviours  are  plotted. 

One  property  of  the  RGA  is  that  the  sum  of  all  the  elements 
of  a  row  or  a  column  is  one  [1 1] .  In  the  case  of  two  dimension 
RGAs,  knowing  one  of  the  four  elements  of  the  matrix,  the 
others  are  quickly  known.  Because  of  that,  only  the  (1,  1) 
element  is  plotted  (Figs.  1  and  2).  As  has  been  explained, 
diagonal  RGA  (RGA(1, 1)  =  1)  are  preferred.  CN  and  MRI  do 
not  depend  on  the  control  structure  pairing  and  in  Figs.  3-6, 
the  control  structure  names  are  indicated  regardless  of  the 
order. 


Fig.  3.  MRI  at  operating  point  one. 
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Fig.  4.  MRI  at  operating  point  four. 


Fig.  5.  CN  at  operating  point  one. 


5.1.  Analysis  at  OP  1 

(1)  In  Fig.  1,  only  23  of  the  36  control  structures  are  rep¬ 
resented.  Absent  structures  are  the  ones  with  the  worst 
behaviour. 

(2)  All  control  structures  with  the  variable  yo2,ca,in  in  the 
set  of  manipulated  variables  have  RGA(1,  1)  very  close 
to  one,  indicating  a  decoupled  multivariable  system.  The 
pairing  in  all  cases  is  such  that  yo2,ca,in  is  in  charge  of 
vst  control.  This  superior  behaviour  is  independent  of  the 
frequency. 


(3)  Three  control  structures  with  FFan,in  controlling  A p 
also  have  a  RGA(1,  1)  very  close  to  one.  They  are 

(0ca,in  Wan,in)?  (^sm,out  Wan,inX  &nd  (Tcm  Wan,in)- 

The  next  group  of  structures  with  well  behaved  RGA 

(1,  1)  are  (0ca,in  —  -^H2,an,in)?  (^sm,out  —  ^H2,an,in)?  and 
(Tcm  —  -^H2,an,in)- 

(4)  In  Fig.  3,  the  best  MRI  behaviours  are  plotted.  ( vcm  — 
Wan, in)  and  (fern  “  *H2,an,in)  have  the  largest  MRI  of 
all  the  control  structures  at  steady  state  and  (yo2,ca,in  — 
Wan 

,in  ),  ( yo2  ,ca,in  -*Tl2,an,in)>  (};02,ca,in  ^rm,outX  and 
(jo2,ca,in  —  Fcm)  have  the  largest  MRI  at  high  frequen¬ 
cies.  It  is  remarkable  that  again,  structures  with  the  vari¬ 
able  yo2,ca,in  are  some  of  the  best  ones  (except  for  the 
(};02,ca,in  0ca,in))-  (Fern  ^rm,out)  has  also  a  good  be¬ 

haviour. 

(5)  With  respect  to  the  CN,  the  best  control  structure  is 
(yo2,ca,in  —  ^rm,outX  followed  by  other  structures  with 
the  variable  yo2,ca,in-  (.y02,ca,in  —  ^ca,out)>  (};02,ca,in  — 
^ca,out)?  and  (yo2,ca,in  —  ^sm,out)-  If  3’o2,ca,in  is  not  con¬ 
sidered,  the  best  control  structures  with  respect  to  the 

CN  are  (.^fsm,out  ^ca  ,out)?  (Fcm  Wan, in),  and  (v  cm 
VH2,an,in)-  The  next  group  of  structures  with  well  behaved 

CN  are  (A"sm,out  ^rm,out)  and  (Fcm  ^frm,out)- 

(6)  Results  show  different  behaviour  of  control  structures 
with  different  location  of  fluidynamic  restriction  CKan,out> 
^ca,out,  ^sm,out,  or  ^rm,out).  This  indicates  that  controlling 
the  system  with  one  Knozz\Q  or  another  has  not  the  same 
effects.  However,  Kc^out  and  Kr m,0ut  will  produce  similar 
behaviour  because  there  is  only  a  small  volume  between 
these  two  fluid  ways. 

Considering  the  results  of  the  three  indexes,  the  following 
general  conclusions  for  the  controllability  around  OP1  are: 

(1)  Control  structures  with  yo2,ca,in  present  the  most 
favourable  indexes:  RGA  closest  to  unity  matrix,  highest 
MRI,  and  smallest  CN. 

(2)  If  yo2,ca,in  is  not  considered,  the  following  best  con- 
trol  structures  are  (ucm  -  Wanjn)  and  (ucm  -  *H2,an,in)- 
These  structures  have  very  similar  behaviour  with  re¬ 
spect  to  the  controllability  indexes.  However,  VH2,an,in 
should  not  be  used  as  manipulated  variable  because  of 
the  delicate  equilibrium  of  the  water  management  sys¬ 
tem. 

(3)  other  structures  with  good  indexes  are  ( vcm  —  ATrrn,out)> 
followed  by  ( vcm  —  /Tan,out)-  The  last  one  has  very  good 
properties  at  low  frequencies,  but  not  at  high  frequencies. 

5.2.  Analysis  at  OP2  and  OP 3 

Analysing  the  MRI  and  CN  in  the  frequency  domain  at 
OP2  and  OP3,  some  changes  are  observed  with  respect  to  the 
analysis  at  OP1 .  Therefore,  the  controllability  indexes  depend 
on  the  linearisation  point.  However,  the  structures  presenting 
the  higher  MRI  are  the  same  as  at  OP1  ((vcm  —  Wan,  in)  > 

(Fcm  —  -^H2,an,in)>  (};02  ,ca,in  Wan,  in  ),  ( yo2  ,ca,in  -Ht2,an,in)> 
(};02,ca,in  —  ^rm,outX  (};02,ca,in  —  FCm)X  2nd  the  MRI 


Fig.  6.  CN  at  operating  point  four. 
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values  do  not  differ  more  than  10%.  The  structures  pre¬ 
senting  smaller  CN  are  also  the  same  as  at  OP1.  Only 
(^sm,out  -  ^ca,out)  has  a  CN  quite  different  in  OP1,  OP2, 
and  OP3.  Therefore,  moving  from  OP1  to  OP2  or  OP3, 
the  sensitivity  of  the  system  is  almost  unchanged  for  these 
preferred  structures.  Finally,  RGA  analysis  indicates  that 
some  input  pairs  would  require  different  pairings  at  OP2 
and  OP3.  This  is  the  case  of  the  pair  formed  by  Kc a?out  and 
fFan?in.  Evidently,  these  structures  could  not  work  properly 
around  both  operating  points.  However,  the  best  structures 
according  to  RGA  are  the  same  at  OP1,  OP2,  and  OP3.  As 
conclusion,  when  the  system  is  controlled  at  constant  ust 
and  A p,  the  preferred  control  structures  will  not  change  due 
to  a  variation  in  the  operating  conditions  in  a  wide  range  of 
operation  and  the  same  tuning  parameters  are  expected  to 
be  adequate  for  operation  around  OP1,  OP2,  and  OP3. 

5.3.  Analysis  at  OP4 

Analysing  the  MRI,  CN  and  RGA  in  the  frequency  range 
between  10-2  and  102rads_1  for  the  system  linearised  at 
OP4  (see  Figs.  2,  4  and  6),  the  following  can  be  stressed: 

(1)  The  three  general  considerations  done  for  the  optimal 
operating  point  stand. 

(2)  Comparing  results  for  OP1  and  OP4,  the  most  interest¬ 
ing  Structures  ((ncm  Wan, in)?  (l?cm  -*n2,an,in)?  (^cm 
^an,out),  and  (vcm  -  Krm,out))  have  more  than  two  times 
better  indexes  at  all  frequencies  for  OP4.  This  result  is  in¬ 
teresting  because  it  indicates  a  possible  trade  off  between 
optimal  operation  (minimal  consumption)  and  controlla¬ 
bility. 

6.  Simulations 

Results  of  Section  5  are  based  on  a  linear  model.  The  non- 
linearities  of  the  system  are  not  considered.  In  order  to  study 
the  nonlinearities  influence,  different  simulations  are  done 
with  the  non  linear  model.  PI  controllers  are  implemented 
in  the  control  loops.  In  Fig.  7,  the  response  of  the  system  to 
changes  in  7st  is  shown.  The  7st  profile  can  be  seen  in  Fig.  7h. 

According  to  Section  5  results,  the  more  interesting 
structures  are  some  of  the  structures  with  yo2,ca,in>  as 
well  as  (vcm  —  Wanjn),  (ncm  —  VH2,an,in)>  (^cm  —  ^rm,out)  and 
(vcm  —  TGn  out)-  From  the  most  interesting  structures  with 
};o2,ca,in  Cyo2,ca,in  -  Wan, in)  has  been  chosen  for  simulation. 
On  the  other  hand,  structures  implying  variables  related  to  the 
system  humidity  are  rejected  because  their  appropriate  values 
have  to  be  maintained  in  order  to  keep  the  delicate  water  equi¬ 
librium.  In  consequence,  the  structures  chosen  for  the  sim¬ 
ulation  comparison  are  (yo2,ca,in  -  Wan?in),  ( vcm  -  Wan?in), 

(^cm  ^rm,outX  and  (ncm  Kan, out)- 

Results  of  Section  5  indicate  that  (yo2,ca,in  —  Wan  in) 
and  ( vcm  —  Wan?in)  have  almost  constant  controllability  in¬ 
dexes  at  OP1,  OP2,  and  OP3.  It  has  also  been  seen  that 


(};o2,ca,in  —  Wan  in)  is  less  affected  by  the  directionality  of  the 
system  and  has  less  influence  between  loops.  As  expected, 
simulations  show  that,  with  unique  tuning  parameters,  both 
control  structures  are  able  to  drive  the  system  from  one  point 
to  another  with  acceptable  performances.  Referring  to  the 
operating  range  limits,  steps  of  ±50%  in  the  current  load 
(not  shown)  are  affordable  for  both  structures  with  smooth 
tunings  (if  the  yo2,ca,in  loop  is  tuned  for  an  aggressive  re¬ 
sponse,  yo2,ca,in  can  arrive  to  its  physical  limits).  However, 
simulations  have  shown  some  differences  between  these  two 
control  structures,  which  have  been  tuned  to  have  similar 
output  peaks.  Since  the  compressor  dynamics  only  affects 
(vcm  —  Wan?in),  this  structure  is  a  bit  slower,  as  can  be  ap¬ 
preciated  in  Fig.  7a  and  b.  The  total  pressure  of  the  cath¬ 
ode  is  much  better  maintained  by  (yo2,ca,in  —  Wan?in)  and 
this  provokes  smaller  changes  in  Wan?in,  as  can  be  seen  in 
Fig.  7e.  Another  advantage  of  (yo2,ca,in  —  Wan?in)  that  can 
be  seen  in  Fig.  7f,  is  that  this  structure  can  maintain  Xo2 
almost  constant.  In  Fig.  7 g,  the  net  power  is  shown.  With 
(vcm  —  Wan?in)  structure,  the  power  consumed  by  the  com¬ 
pressor  is  higher  when  the  7st  is  higher  and  lower  when  the  7st 
is  lower  and  therefore,  the  net  power  increments  and  decre¬ 
ments  are  always  diminished  resulting  in  a  more  constant 
value.  Since  (yo2,ca,in  -  Wan,in)  does  not  modify  vcm ,  the 
power  consumed  by  the  compressor  is  almost  constant  and 
the  net  power  more  variant  with  the  7st  changes. 

When  analysing  the  control  structures  ( vcm  —  Km,out)  and 
(vcm  —  Krm ?out),  very  different  behaviours  are  found  close 
to  the  nominal  operating  point  and  far  from  this  operating 
point.  These  structures  are  not  able  to  control  the  system  if 
the  load  current  has  large  variations  because  of  the  hydrogen 
exhaustion.  In  Fig.  7a  it  is  seen  that  ( vcm  —  7frm?0 ut)  can 
afford  the  first  current  step  but  is  not  able  to  control  the  sys¬ 
tem  for  the  second  current  step.  ( vcm  —  7fan?out)  has  a  wider 
range  of  operability  than  ( vcm  —  7frm?ou t)  but  its  operability 
range  is  still  much  smaller  than  the  operability  range  of 
(x>2,ca,in  -  wan, in)  and  ( vcm  -  Wm,m)-  Therefore,  these  two 
structures  would  work  properly  close  to  the  nominal  operat¬ 
ing  point  but  would  not  be  able  to  drive  the  system  for  large 
excursions. 

Simulation  results  show  peaks  in  ko2  (Fig.  If).  These 
peaks  are  inadequate  and  can  provoke  oxygen  starvation  and 
the  rapid  degradation  of  the  membrane.  Through  simulations 
which  are  not  shown,  Xo2  behaviour  has  been  compared  in 
two  different  cases:  the  control  of  nst  by  vcm  and  the  control 
of  \o2  by  vcm.  Results  indicate  that  similar  Xo2  peaks  are  ob¬ 
tained.  Tuning  parameters  hardly  influence  the  value  of  this 
peak.  Therefore,  the  control  of  Xo2  directly  cannot  solve  the 
problem  of  the  quick  peak  in  Xo2 .  In  addition,  the  measure 
of  Xo2  would  be  difficult. 

7.  Discussion 

Results  of  linear  analysis  indicate  that  control  structures 
using  yo2,ca,in  would  be  very  adequate,  and  simulations 
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with  structure  (yo2,ca,in  —  ^an,in)  confirm  its  appropriate 
behaviour.  However,  manipulation  of  yo2,ca,in  would 
imply  an  added  complexity  (N2  and  O2  sources  would 
be  required).  This  added  complexity  will  compensate  or 


not  the  controllability  improvement  in  dependence  of  the 
application. 

As  has  been  seen,  when  the  fuel  cell  is  operated  close  to 
a  nominal  operating  point,  control  structures  that  use  ^an,out 
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Fig.  7.  Simulation  results  for  7st  variation  shown  in  h. 
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are  possible.  Using  this  control  variable,  hydrogen  output 
flow  varies  with  7st.  On  the  contrary,  if  Wa n3n  is  used,  hy¬ 
drogen  output  flow  is  better  maintained  at  the  nominal  value, 
what  is  important  for  a  regular  behaviour  at  the  end  of  the 
distribution  channels.  ( vcm  —  Wa n,in)  is  clearly  superior  than 
(vcm  —  Kan  ouf),  but  this  second  structure  can  still  be  useful 
in  case  of  failure  in  the  ( vcm  —  Wan,in)  control  loop,  for  a  cor¬ 
rect  shut  down,  for  example.  Close  to  a  nominal  operation, 
also  ( vcm  —  ^rm,out)  is  a  valid  structure.  However,  the  limita¬ 
tions  in  this  case  are  still  more  restrictive.  Being  both  control 
variables  in  the  cathode  side,  conflicts  between  loops  appear 
easily. 

A  rapid  control  of  ust  appears  to  be  contradictory  with 
slow  variations  of  po2  because  of  the  close  relation  between 
these  two  variables.  Therefore,  it  seems  that  in  order  to  avoid 
po2  peaks  and  oxygen  starvation,  vsi  settling  time  must  be 
relaxed.  This  limit  of  the  fuel  cell  dynamic  behaviour  will 
influence  the  design  of  the  conditioning  system  and  its  control 
objectives. 


8.  Conclusions 

Through  linear  analysis  tools,  the  controllability  of 
different  control  structures  for  polymeric  fuel  cells  has  been 
compared.  Results  of  the  linear  analysis  indicate  that  different 
structures  are  suitable,  what  is  confirmed  through  simulation. 
Several  control  structures  are  adequate  if  there  is  a  tight  oper¬ 
ation  around  a  nominal  operating  point,  which  is  possible  if 
the  fuel  cell  is  provided  with  a  good  electric  power  condition¬ 
ing  system.  However,  as  has  been  seen  through  simulations, 
not  all  these  structures  are  valid  when  the  operating  condi¬ 
tions  move  far  from  the  nominal  conditions.  ( vcm  —  fkan  in) 
and  (yo2,ca,in  —  Wan, in)  control  structures  have  the  largest 
operating  ranges.  Of  these  two  structures,  (yo2,ca,in  —  Wan?in) 
is  less  affected  by  the  directionality  of  the  system  and  has  less 
influence  between  loops.  The  use  of  yo2,ca,in  for  the  control 
of  ust  has  other  important  advantages  that  could  justify 
the  additional  complexity  required  to  change  the  oxidant 
composition.  Simulations  have  also  shown  that,  although  a 
limitation  on  the  vst  settling  time  is  necessary,  fuel  cells  can 
be  controlled  by  linear  decentralised  feedback  controllers 
with  PI  in  each  loop.  Comparing  the  controllability  indexes 
of  a  polymeric  fuel  cell  at  different  operating  conditions  it 
has  been  found  that  a  trade  off  between  control  performance 
and  hydrogen  consumption  is  possible.  The  controllability 
of  the  system  can  be  improved  changing  the  operating 
point. 
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Appendix  A 


The  equations  for  the  fuel  cell  voltage  are: 


E  =  1.229  -  8.5  x  l(T4(Tfc  -  298.15)  +  4.3085 


x  10_57fc 


1 

ln(/>H2)+  2ln(/>o2) 


fact  =  VQ  +  fa(l  -  ecu) 


i>o  =  0.279  -  8.5  X  10_4(7fc  -  298.15)  +  4.3085  x  10“+ 


X 


In  (  ^ca  ^sat^\  .  1  i  { 0.1173 (Pea  j^sat)^ 

n  1.01325  )  +  2  {  1.01325  J  J 


fa  = 


(-1.618  X  10“+  +  1.618  X  10“+ 

V0.1173 

+(1.8  x  10“+  -  0.166)  (  P°2 


Po2  x2 


+  P  sat 


0.1173 


+  Psat 


+(-5.8  x  10_47fc  +  0.5736) 


^ohm  —  i  x  ^ohm 


^ohm  — 


'm 


a, 


m 


=  (fciAm  -  bn)  exp  b2 


1 


1 


303  7fc 


C3 


^conc  —  M  c2  . 

I 


max 


with  the  following  values  of  the  experimental  parameters: 


'it  f)(h 


0.1173 


+  Psat  <  2, 


C  2  =  { 


(7.16  x  10_47fc  -  0.622)  +  pj) 

+(—1.45  x  10~37fc  +  1.68)  else, 

(8.66  x  10_57fc  -  0.068)  +  pj) 


+(-1.6  x  10-4Te  +  0.54) 


a  =  10 


b  ii  =  0.005139 


b  12  =  0.00326 


b  2  =  350 

6nax  =  2. 

c3  =  2 


102 
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where  tm  is  the  membrane  thickness,  and  Am  is  the  membrane 
water  content. 


Appendix  B 


Water  transport  across  the  membrane  Wv,membr  is 

tt j  Ti/r  a  i  ^d^st  ^  ^v,ca  —  cv,an 

Wv?membr  —  MyA{cfl  (  —  7)w 


F 


t 


m 


d[  =  ,  i  =  [anode,  cathode] 

Ps2Lt,l 


Am  — 


ifO  <  d{  <  1,  0.043  +  17.81^1  —  39.85 a?  +  36a? 

if  1  <  di  <3,  14  +  1.4(ai  —  1) 


nd  =  0. 00291^  +  0.05Am  -  3.4  x  10“19 
Dw  =  D, exp  (2416(^-2-)) 

'  if  Am  <  2,  10“6 

if 2  <  Xm  <  3,  10-6(l+2(V-2)) 

Oa  =  < 

if  3  <  Am  <  4.5,  10-6(3  —  1.67(Xm  —  3)) 

if  A.m  <  4.5,  10~6  x  1.25 

where  cv,c a  and  cv,a n  are  the  water  concentrations  in  cathode 
and  anode  sides,  and  is  the  molar  fraction  of  water. 
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